Parametrized canonical transformation for the Hubbard-model at arbitrary 

interaction strength 
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The t — J and Heisenberg models are truncated expansions of a canonically transformed Hubbard 
model coinciding with it at U — > oo. We show that a modified canonical transformation applied to 
the Hubbard model leads to alternative models of similar form, but whose convergence properties 
with respect to the expansion are more favourable, resulting in a good description of the half-filled 
ground state even at < U < 1. We investigate the transformed Hamiltonian and observables for 
metallic and insulating variational wave-functions. 

PACS numbers: 71.10.Fd,71.30+h 



The Hubbard model [H, 0, H, 0] and its descendants 
have contributed greatly to our understanding of strongly 
correlated systems (51 H @] and in particular the metal- 
insulator transition^] (MIT) exhibited by these systems. 
Early attempts d @] to explain the MIT were based on 
the use of a projected wavefunction due to Gutzwillcr 
(GW). An approximate variational calculation based on 
the Gutzwiller approximation(GA) [J] for the GW in the 
general case predicts a MIT |8j between a paramagnetic 
metal and an insulator (Brinkman-Rice transition). The 
order parameter for the Brinkman-Rice transition is the 
fraction of doubly occupied sites which goes to zero at the 
critical U = U c . A shortcoming of the GA is that second- 
order hopping processes are not included, i.e. double oc- 
cupations that arise as a result of second-order hoppings 
(which give rise to anti-ferromagnetic (AFM) coupling) 
are entirely absent. Thus the number of double occupa- 
tions is not a valid order parameter for the actual MIT. 
In one dimension, the exact solution for the Hubbard 
model fj| indicates insulating behavior for all finite val- 
ues of the interaction, whereas the exact solution for the 
GW [l(| for the same system is always metallic. 

The importance of higher-order hopping processes is 
made obvious by a canonical transformation of the Hub- 
bard model that eliminates those first-order hopping pro- 
cesses which increase (decrease) the number of doubly oc- 
cupied sites (iJ+(ff t ")) p], El El El EEL El- Expan- 
sion and truncation of the transformed Hamiltonian leads 
to the well-known t — J and spin— ^ anti- ferromagnetic 
Heisenberg models, which coincide with the Hubbard 
model in the strong- coupling limit, in which it leads to 
anti-fcrromagnetism. [lot Il7| 

The effective Hamiltonians derived from the Hubbard 
model have other applications as well. In the resonating 
valence bond (RVB) method [H, EH, H3] the expecta- 
tion value of the t — J Hamiltonian is evaluated over a 
fully Gutzwiller projected wavefunction 0,13] . The RVB 
wavefunction has recently been applied to the problem 
of high temperature superconductivity, and many exper- 
imentally observed features of the relevant materials have 
been reproduced. @, El, SI] 

In the present study the unitary operator that trans- 
forms the Hubbard model into the t — J or Heisenberg 



models is parametrized so that the number of double oc- 
cupations as a function of the transformation can be 
minimized. The effect of our procedure is similar to 
that of the original transformation. The difference is 
that -ff t + and Hj~ are not cancelled from the Hamilto- 
nian as in the standard case, but instead constrained 
so that their expectation values are zero. In contrast, 
the t — J and spin-^ Heisenberg models will in general 
give finite expectation values for _ff t + and H^~ . In our 
approach first-order double occupations are eliminated 
at the wavefunction level, as opposed to the operator 
(Hamiltonian) level. The optimized transformation can 
be applied at any value of the interaction and not only in 
the strongly interacting limit. We diagonalize the trans- 
formed Hamiltonians for systems of up to 12 lattice sites, 
and it is shown that the optimized expansion converges 
much faster than the standard one. Convergence is also 
demonstrated for U < 1. 

We also investigate the behavior of the optimally trans- 
formed double occupation operator using two different 
variational wavefunctions the Gutzwiller 3 (GW) and 
Baeriswyl [23] (BW) wavefunctions and compare them 
to the exact result. 

The Hubbard model Hamiltonian can be written as 



Hu 

UD 



(1) 



where D = X)i n »T n il an< ^ where the operator c ia .(ci a ) 
creates(destroys) a particle at site i with spin er, and ni a 
is the density operator at site i for particles of spin a. 
In deriving the canonically transformed Hamiltonian it 
is helpful to break up the kinetic energy operator into 
terms consisting of different types of hoppings [7[ : 



H t = H+ + 



H 1 



(2) 
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where 

H+ = -t ^ m-achcjvil - Uj-a) (3) 

-i ^ (1 - n i _ CT )cj (T c J(r (l - rij-v) 

H^(H^r) include only hopping processes which in- 
crease(decrease) the number of double occupations, and 
H® includes only those which leave the number of double 
occupations unchanged. The Hermitian operator defined 
as 

S=-^{H+-Ht) (4) 
is useful in defining the transformation 

H S = e lS He- lS = H + i[S,H] + l - [S, [S, H}} + ... (5) 

The series can be viewed as a power series in jj. It can 
be shown that 

i[S,H u ] = -(H++H t ~), (6) 

and thus, up to first order, hoppings that change the 
number of double occupations are cancelled from the 
transformed Hamiltonian (Eq. 1(5))). The t — J and 
Heisenberg models, which are used as effective models in 
the large U limit, can be derived by explicitly evaluating 
the terms of Eq. ([5J up to second order in t/U, 

H S « fl? + H V + J£ { . j} (Si • S, - ^) (7) 
+3-site terms 

where J = 4t 2 /U. 

We now consider a similar transformed Hamiltonian 
derived using the modified operator e % which leads to 

H aS = e iaS He- iaS (8) 

■2 2 

= H + ia[S,H] + l -^-[S,[S,H]} + ..., 

where a is a parameter to be determined. If for a partic- 
ular state the transformed number of double occupations 

(y\D aS \^) = (y\e iaS De- iaS \V), (9) 
is minimized as a function of a, then it holds that 

(y\e laS [S,D}e- iaS \^} = 0, (10) 
which with Eq. ([6]) is equivalent to 

(ty\e iaS {H+ + Hf)e- iaS \V) = 0. (11) 
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FIG. 1: Optimal a as a function of U for systems with differ- 
ent sizes. 



Hamiltonian 


U 


2nd order 


4th order 


6th order 


Exact 


e lS He- lS 


0.5 


-395.505 


-4613.096 


-35947.499 


-7.275 




1.0 


-99.211 


-270.19 


-495.743 


-6.601 




2.0 


-24.713 


-10.3987 


-7.6983 


-5.409 




5.0 


-4.557 


-2.974 


-3.092 


-3.088 




10.0 


-1.824 


-1.661 


-1.664 


-1.664 


e iaS He~ iaS 


0.5 


-11.084 


-6.695 


-7.328 


-7.275 




1.0 


-9.850 


-6.159 


-6.634 


-6.601 




2.0 


-7.742 


-5.158 


-5.421 


-5.409 




5.0 


-3.819 


-3.047 


-3.088 


-3.088 




10.0 


-1.792 


-1.662 


-1.664 


-1.664 



TABLE I: Comparison of ground state energies calculated for 
a lattice composed of six sites. The upper(lower) half shows 
results for the transformed Hamiltonian with a = 1 (optimized 
a). The rightmost column shows the exact results. The ex- 
pansion is in the parameter a. 



Thus, double occupations up to first-order can be ex- 
cluded via a transformation that sets the expectation 
value of the sum of the operators H t + + to zero. The 
main difference between the Hamiltonians in Eq. (J5J> and 
Eq. (|8|) is that in the latter the expectation value of the 
sum of the operators that change the number of double 
occupations is zero, as opposed to being cancelled by an- 
other term equal but opposite in sign at the operator 
level. 

If $ is the ground state of the Hubbard Hamiltonian, 
then 

(&\H\$) = <$«s|tf Q s|*as}, (12) 

where the transformed wavefunction |$ Q s) = e* |$) is 
the ground state of the transformed Hamiltonian H a s- 
While the optimization procedure can be carried out on 
any state, in the rest of this work we deal exclusively with 
the ground state at half filling. 

The analog derivation that leads to the t — J model 



3 



applied to Eq. ([8]) results in 

H aS a H? + H U + J aS (Si ■ Sj - (13) 

+3-site terms, 

where J Q s denotes a modified coupling constant satisfy- 
ing 

J aS = {2a - a 2 ) J. (14) 

The first-order term in a originates from the transformed 
H+ and Ht. 

The size of the parameter a determines the con- 
vergence of the expansion (Eq. ©). In Fig. Q] 
the results of power-method type calculations [26j are 
shown for systems of various sizes at half-filling. Anti- 
periodic (periodic) boundary conditions were ap plie d for 
system sizes with odd(even) multiples of two [27|, l28j . 
The parameter a which minimizes Eq. ^ ( and satis- 
fies Eqs. (|10p and (JTTJ)) and is closest to the origin is 
calculated as a function of the interaction parameter U. 
We find that convergence is achieved for all U consid- 
ered. As expected, Hs is recovered for large U. The 
size-dependence of a is negligible. Interestingly, as U 
approaches zero a/U converges to « 0.3, wheras in the 
standard case 1/U diverges. 



to the exact result in all cases, and the convergence is also 
better when the expansion of the Hamiltonian is carried 
out to higher orders. The advantage is more pronounced 
at lower values of U, in particular our transformation is 
even applicable for U < 1 where the standard expansion 
fails due to slow convergence. The second order results 
with optimal a (similar to the t — J model) are in consid- 
erably better agreement with the exact results than the 
standard (a = 1) second order ones, therefore the t — J 
model is, in this sense, applicable even at U < 1, but 
with a modified coupling. 
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FIG. 2: Ground state expectation value of the transformed 
interaction U (D a s)/t for the standard expansion and the 
optimized one compared to the exact results for a model with 
six sites. The transformed D a s was expanded to second order 
in a in the main plot, sixth order in the inset. 

In Table H] we compare energies calculated using the 
standard transformation (Eq. ©), and those resulting 
from the transformation with optimized a (Eq. ([8]) and 
Fig. [1]) . The optimal value of a was obtained from exact 
diagonalization. In these calculations periodic boundary 
conditions were used. Subsequently, a was used in the 
expansion, Eq. ([8]). In order to investigate the conver- 
gence, the expansion of the Hamiltonian was carried out 
to second, fourth, and sixth orders in a, then diagonal- 
ized. The optimized transformation gives energies closer 



FIG. 3: Ratio SI (defined in Eq. ||15| l1 calculated exactly for 
different system sizes. The inset shows a comparison between 
the exact result and two different variational wavefunctions 
(Baeriswyl (BW) and Gutzwiller (GW)) for the system with 
12 lattice sites. 

In Fig. [2] the expectation value of the transformed in- 
teraction energy is shown. The expansion is carried out 
to second and sixth order (inset) for a = 1 and for opti- 
mized a, i.e. the Haniltonian is calculated up to a given 
order, and diagonalized. The observable is also trans- 
formed and truncated at the given order. Optimized a 
gives quantitative agreement with the exact result even 
at second order (t — J like model) , whereas the standard 
version is not in agreement with the exact results at sec- 
ond order, and even when the expansion is carried out to 
sixth order, agreement is only reached when U is large. 

The t — J type model derived herein is not as easy 
to derive as the standard one. At a particular U the 
normal t — J model can easily be derived to any order. 
Our modified model depends on a parameter, a, which 
is a function of the ground state solution. For a partic- 
ular U one can obtain a by expanding the transformed 
Hamiltonian (Eq. (j5J|), solving for its ground state, and 
varying a to satisfy the condition in Eq. (fTU)) . It also 
appears possible to apply our formalism using the gen- 
eralized version of the canonical transformation of Ref. 
0. 

We have also investigated our scheme for different vari- 
ational wavefunctions. For our studies we have chosen 
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the Baeriswyl and Gutzwiller wavefunctions (BW and 
GW respectively). The properties of these wavefunctions 
are well-known. In particular it has been shown by Millis 
and Coppersmith [2J] that the Drude weight of the GW 
is always finite in the thermodynamic limit, hence the 
GW is metallic. This property can be attributed to the 
lack of explicit phase dependence of the GW. The BW 
has been shown to consist of rotating dipoles formed of 
empty and doubly occupied sites, and to be in general an 
insulating wavefunction (25|. 

In Fig. [3] we present a comparison of the ratio 



n = 



(*|L>|*) 



(15) 



for systems with different sizes calculated exactly. As U 
increases ft decreases sharply. The inset in Fig. [3] shows 
a comparison for the system of size 12 between the exact 
result and two variational wavefunctions BW and GW. 
An interesting feature is that in the large U limit the 
GW tends to a finite value unlike the exact or the BW 
result. These qualitative tendencies persist away from 
half- filling (results not shown). Hence the GW tending 
to a finite limit is not due to metallicity. 
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FIG. 4: Optimal a as a function of U for the Gutzwiller 
wavefunction with 12 sites for different fillings. 
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In Fig. [4] we show the optimum a at three different 
fillings for the GW. At half-filling the behaviour is quali- 
tatively different from the other fillings investigated, and 
different from the behaviour found for the exact case (Fig. 
[IJ. At large U a appears to be bounded below for half- 
filling, where GW is expected to be in error, since it is a 
metallic wavefunction. Away from half-filling the a ob- 
tained from GW is monotonically increasing. We have 
also investigated the BW and found the qualitative ten- 
dencies (monotonic increase, upper bound of a = 1) to 
be the same as for the exact calculation. 



In conclusion we have shown that the standard canon- 
ical transformation which when applied to the Hubbard 
model gives the t — J model at large interaction strength 
can be optimized to give a i — J like model applicable for 
the whole range of the interaction strength. In particu- 
lar convergence of the expanded Hamiltonian is achieved 
for interaction strength close to zero, where the standard 
transformation leads to slow convergence. 



Beneficial discussions with E. Arrigoni and W. von der 
Linden are gratefully acknowledged. 
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